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Abstract 

The prediction of the effective elastic properties of polymer bonded explosives using direct numerical simulations is 
computationally expensive because of the high volume fraction of particles in these particulate composites (~0.90) 
and the strong modulus contrast between the particles and the binder (~20,000). The generalized method of cells 
(GMC) is an alternative to direct numerical simulations for the determination of effective elastic properties of compos- 
ites. GMC has been shown to be more computationally efficient than finite element analysis based approaches for a 
range of composites. In this investigation, the applicability of GMC to the determination of effective elastic properties 
of polymer bonded explosives is explored. GMC is shown to generate excellent estimates of effective moduU for com- 
posites containing square arrays of disks at volume fractions less than 0.60 and a modulus contrast of approximately 
100. However, for high volume fraction and strong modulus contrast polymer bonded explosives such as PBX 9501, 
the elastic properties predicted by GMC are found to be considerably lower than finite element based estimates and 
experimental data. Simulations of model microstructures are performed to show that normal stiffnesses are underes- 
timated by GMC when stress-bridging due to contact between particles is dominant. Additionally, the computational 
efficiency of GMC decreases rapidly with an increase in the number of subcells used to discretize a representative 
volume element. The results presented in this work suggest that GMC may not be suitable for calculating the effective 
elastic properties of high volume fraction and strong modulus contrast particulate composites. Finally, a real-space 
renormalization group approach called the recursive cell method (RCM) is explored as an alternative to GMC and 
shown to provide improved estimates of the effective properties of models of polymer bonded explosives. 

Keywords : Effective Properties, Particulate Composites, High Volume Fraction, Strong Modulus Contrast, Stress- 
Bridging, Method of Cells, Real-Space Renormalization Group 

1 Introduction 

The generaUzed method of cells (GMC) (Aboudi 1996, Paley & Aboudi 1992) is a semi-analytical method of deter- 
mining the effective properties of composites. In this method, a representative volume element (RVE) of the composite 
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under consideration is discretized into a regular grid of subcells. Equilibrium and compatibility are satisfied on an av- 
erage basis across subcells using integrals over subcell boundaries. GMC generates a matrix of algebraic expressions 
containing information about subcell material properties. The effective stiffness of the composite can be obtained by 
inverting this matrix. 

One advantage of GMC over other numerical techniques is that the full set of effective elastic properties of a com- 
posite can be calculated in one step instead of solving a number of boundary value problems with different boundary 
conditions. GMC has also been found to be more computationally efficient that finite element calculations for fiber 
reinforced composites (Aboudi 1996, Wilt 1995), since far fewer GMC subcells than finite elements are necessary to 
obtain the same degree of accuracy. The problem of discretization is also minimized since a regular rectangular grid 
is used in GMC. 

The generalized method of cells is discussed briefly in this work. Effective stiffnesses predicted by this method are 
compared with accurate numerical predictions for square arrays of disks (Greengard & Helsing 1998). The method is 
then applied to two-dimensional models of a general polymer bonded explosive and to the microstructure of PBX 950 1 
using a two-step procedure similar to that of Low et al. (Low et al. 1994). GMC estimates of elastic properties are 
compared with predictions from detailed finite element calculations. The performance of GMC is explored for several 
microstructures with contacting particles and some shortcomings of the method are identified. Finally, the recursive 
cell method (RCM) is explored as an alternative to direct GMC calculations in the prediction of effective properties of 
polymer bonded explosives. 



2 The generalized method of cells 

Figure [T] shows a schematic of the RVE, the subcells and the notation (Aboudi 1991) used in GMC. In the figure, 
(Xi, X3) is the global coordinate system of the RVE and (x'f\x2^\x'"^'') is the coordinate system local to a subcell 
denoted by (a/Sj). It is assumed that the displacement function varies linearly within a subcell (a/S^) and can 

be written in the form 
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where i represents the coordinate direction and takes the values 1, 2 or 3; w\" ' is the mean displacement at the 
center of the subcell (a/37); and ^^°'^'^\ and ^^"^'•'^ are constants local to the subcell that represent gradients 

of displacement across the subcell. The strain-displacement relations for the subcell are given by 

e(f^Ul(a,uf^^'+a,-nf^^') (2) 

where di = d/dx^i \ 82 = d/dx'f \ and 83 = d/dx^\ Since polymer bonded explosives are isotropic particulate 
composites, the following brief description of GMC assumes that the RVE is cubic and all subceUs are of equal size. If 
each subcell (a^-y) has the same dimensions (2/i, 2h, 2h) then the average strain in the subcell is defined as a volume 
average of the strain field over the subcell as 

The average strain in the subcell can be obtained in terms of the displacement field variables. It is assumed that there 
is continuity of traction at the interface of two subcells. The displacements and tractions are assumed to be periodic 
at the boundaries the RVE. Applying the displacement continuity equations on an average basis over the interfaces 
between subcells, the average strain in the RVE can be expressed in terms of the subcell strains. The average subcell 
stresses can be obtained from the subcell strains using the traction continuity condition and the stress-strain relations 
of the materials in the subcells. A relationship between the subcell stresses and the average strains in the RVE is thus 
obtained. 

For orthotropic, transversely isotropic or isotropic materials, the approach discussed above leads to the decoupling 
of the normal and shear response of the RVE. This decouphng leads to two systems of equations relating the subcell 
stresses and the average strains in the RVE. For the normal components of strain, the system of equations can be 
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where N is the number of subcells per side of the RVE. The corresponding system of equations for the shear compo- 
nents is of the form 
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In equations Q and (|5]l the M matrices contain material comphance terms. The T matrices contain the average 
subcell stresses. The vector H contains the dimensions of the subcells. Thus, a sparse system of equations of size 3-/V^ 
is produced that relates the subcell stresses to the average strains in the RVE. After inverting these equations and with 
some algebraic manipulation, explicit algebraic expressions for the individual terms of the effective stiffness matrix 
can be obtained. These stress-strain equations that relate the average RVE stresses to the average RVE strains are of 
the form 
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where Cf^ are the terms of the effective stiffness matrix. Details of the algebraic expressions for these terms have 
been pubUshed by other researchers (Pindera & Bednarcyk 1997). 

In GMC, the number of equations to be solved equals the number of subcells raised to the d th power, where d is 
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the number of dimensions in the problem. As a result, the computational efficiency of GMC decreases as the number 
of subcells increases. This issue has been partially resolved (Orozco 1997) by identifying the sparsity characteristics 
of the system of equations and by using the Harwell-Boeing suite of sparse solvers. The computational efficiency 
of GMC has been further improved after a reformulation (Pindera & Bednarcyk 1997, Bednarcyk & Pindera 1997) 
that takes advantage of the continuity of tractions across subcells to obtain a system of 0{N'^) equations in three 
dimensions. 

Due to decoupling of the normal and shear response of the RVE, the shear components of the stiffness matrix 
obtained from GMC are the harmonic means of the subcell shear stiffnesses and of the form 

N N N 

l/Cf, = l/N^Y.T.T.y^ef'^- (7) 

a=l /3 = 1 7=1 

Bednarcyk and Arnold (2001) suggest that this lack of coupUng makes for an "ultra-efficient" micromechanics model. 
However, this lack of coupling can lead to gross underestimation of shear moduli for high volume fraction and high 
modulus contrast materials such as polymer bonded explosives. Recently, researchers (Williams & Aboudi 1999, Gan 
et al. 2000) have attempted to solve the problem by using higher order expansions for the displacement and by ex- 
phcitly satisfying both subcell equilibrium and compatibility. However, these approaches decrease the computational 
efficiency of GMC considerably and are not explored in this work. 



3 Validation - square arrays of disks 



In this section, estimates of effective properties from GMC are compared with accurate numerical results for square 
arrays of disks. Square RVEs containing square arrays of disks exhibit square symmetry. The two-dimensional linear 
elastic stress-strain relation for these RVEs can be written as 
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where Keg is the two-dimensional effective bulk modulus, fi^^ is the effective shear modulus when a shear stress is 

(2) 

applied along the diagonals of the RVE, and /i^jj is the effective shear modulus when a shear stress is applied along the 
edges of the RVE. These three effective moduli have been determined accurately, using an integral equation approach, 
by Greengard and Helsing (1998) for square arrays of disks containing a range of disk volume fractions. 

To compare the effective moduli predicted by GMC with those from the integral equation calculations (Greengard 
& Helsing 1998), RVEs containing disk volume fractions from 0. 10 to 0.70 were created. These RVEs were discretized 
into 64x64 equal sized subcells. The effective stiffness matrix of each RVE was calculated using GMC. Finally, the 
two-dimensional effective moduli for each RVE were calculated from the effective stiffness matrix using the relations 

= OMCf, + Cfi) , fi^^J = 0.5(Cif - Cjf ) , = (9) 

Figure |2] shows the moduli predicted by GMC and those from the integral equation method of Greengard and 
Helsing (G&H) for disk volume fractions from 0.10 to 0.70. The material properties of the disks and the binder used 
in the calculations are shown in Table [T] The effective bulk moduli (Xeff) and diagonal shear moduli (/^gjf ) obtained 
from the GMC calculations are within 4% of those obtained by the integral equation method for all volume fractions 
up to 0.60. At a volume fraction of 0.70, the GMC predictions for bulk modulus and diagonal shear modulus are 

(2) 

4% and 11% less, respectively. For the shear modulus ^^ff , the GMC predictions are around 4% to 10% less than 
the estimates of Greengard and Helsing for volume fractions from 0.10 to 0.60. The difference is around 24% for a 
volume fraction of 0.70. 

These results show that GMC estimates are quite accurate for composites containing square arrays of disks with 
volume fractions up to 0.60, confirming results reported elsewhere (Aboudi 1996). In the next section, GMC is used 
to determine the effective properties of models of polymer bonded explosives and the results are compared to detailed 
finite element calculations and experimental data. 
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4 Modeling polymer bonded explosives 

Polymer bonded explosive (PBX) materials typically contain around 90% by volume of particles surrounded by a 
binder. The particles consist of a mixture of coarse and fine grains with the finer grains forming a filler between 
coarser grains. Modeling the microstructure of these materials is difficult due to the complex shapes of HMX particles 
and the large range of particle sizes. Two-dimensional approximations of the microstructure of PBXs based on digital 
images (Benson & Conley 1999) have been used to study some aspects of the micromechanics of PBXs. However, 
such microstructures are difficult to generate and require complex image processing techniques and excellent image 
quality to accurately capture details of the material. A combination of Monte Carlo and molecular dynamics techniques 
have also been used to generate three-dimensional models of PBXs (Baer 2001). Microstructures containing spheres 
and oriented cubes have been generated using these techniques and appear to represent PBX microstructures well. 
However, the generation of microstructures using dynamics -based methods is extremely time consuming when tight 
particle packing is required, as is the case for volume fractions above 0.70. 

Comparisons of finite element predictions with exact relations for the effective properties of composites (Banerjee 
2002) have shown that detailed finite element estimates can be used as a benchmark to check the accuracy of predic- 
tions from GMC. In this investigation, manually generated PBX microstructures containing symmetrically distributed 
circular particles are used initially to compare GMC and finite element predictions. The two-dimensional microstruc- 
tures contain 90% particles by volume and use two particle length scales. Two-dimensional models containing ran- 
domly distributed circular particles that reflect the actual particle size distribution of PBX 9501 are next modeled with 
GMC and the results compared to finite element estimates. 

The material properties used for the particles, the binder, and PBX9501 in these calculations are shown in Table |2] 
These properties correspond to those of HMX (the explosive particles), the binder, and PBX 9501 at 25° C and a strain 
rate of 0.05/s (Wetzel 1999). 
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4.1 Simplified models of PBX materials 

GMC and finite element calculations were performed for the six, manually generated, simplified model microstructures 
of polymer bonded explosives shown in Figure[3] These representative volume elements (RVEs) contain one or a few 
relatively large particles surrounded by smaller particles to reflect common particle size distributions of PBXs. The 
volume fraction of particles in each of these models is around 0.90±0.005. The binder material surrounds all particles 
in the six microstructures. 

For the GMC calculations, a square grid was overlaid on the RVEs to generate subcells. Two different approaches 
were used to assign materials to subcells before the determination of effective properties of the RVE. In the first 
approach, referred to as the "binary subcell approach", a subcell was assigned the material properties of particles if 
more than 50% of the subcell was occupied by particles. Binder properties were assigned otherwise. Figure|4ja) shows 
a schematic of the binary subcell approach. In the second approach, called the "effective subcell approach", a method 
of cells calculation (Aboudi 1991) was used to determine the effective properties of a subcell based on the cumulative 
volume fraction of particles in the subcell (Banerjee et al. 2000). Figure|4fb) shows a schematic of the effective subcell 
approach. After the subcells were assigned material properties, the GMC technique was used to compute the effective 
properties of the RVE. 

Note that the particles are not resolved well when materials are assigned to subcells in this manner if the number 
of subcells is small. However, the large size of the matrix to be inverted in GMC limits the number of subcells 
that can be used to discretize the RVE. If the binary subcell approach is used to assign subcell materials, contacting 
particles are created where there are none in the actual microstructure, leading to the prediction of higher than actual 
stiffness values. The effective subcell approach improves upon the binary subcell approach by "smearing"the material 
properties at the boundaries of particles and thus reducing the particle contact artifacts caused by discretization errors. 

For validating the GMC results, detailed finite element (FEM) calculations were performed using six-noded trian- 
gular elements to accurately model the geometry of the particles. Around 65,000 nodes were used to discretize each of 
the models. The volume average stress and strain in each RVE was determined for applied normal and shear displace- 
ments. Periodicity was enforced through displacement boundary conditions. Since these finite element calculations 
serve to validate the GMC calculations, further mesh refinement was explored and the results were found to converge 
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those from the 65,000 node finite element models. 

Table [3] lists the effective stiffnesses of the six RVEs shown in Figure [3] from GMC and FEM calculations. On 
average, the GMC calculations using the binary subcell approach predict values of Cfl that are around 2.5 times the 
FEM based values. The values of Cff from the effective subcell approach based GMC calculations are closer to the 
FEM estimates than the binary subcell based GMC estimates. The GMC and FEM estimates of Cjf are quite close. 
The values of C^f from GMC are only 10% of the FEM values. The low values of C^f are obtained because GMC 
predict effective shear stiffnesses that are harmonic means of subcell shear stiffnesses. Models 3 and 4 (Figure[3]l pro- 
duce agreement in Cl^ and Cjf between the binary subcell approach based GMC calculations and FEM (within 5%) 
whereas the other four models produce considerable differences in predictions. These large differences are produced 
by discretization errors introduced in the GMC approach that lead to continuous stress-bridging paths across the RVEs 
and hence to increased stiffness. 

However, if the predicted effective stiffnesses shown in Table |3] are compared with the experimental effective 
stiffness of PBX 9501 (shown in Table|2]i, it can be observed that the models predict values of Cff and Cff that are 
around 10% of the experimental values. Hence, these simplified models are not appropriate for the modeling of PBX 
9501. The next section explores models based on the actual particle size distribution of PBX 9501. 

4.2 Models of PBX 9501 

Coarse and fine particles of HMX are blended in a ratio of 1 : 3 by weight and compacted in the process of manufacturing 
PBX 9501. Figure |5l a) shows four RVEs of PBX 9501 based on the particle size distribution of the dry blend of 
HMX (Wetzel 1999) prior to compaction. Figure |5lb) shows four RVEs based on the particle size distribution of 
pressed PBX 9501 (Skidmore et al. 1998). The larger particles are broken up in the pressing process leading to a 
larger proportion of smaller particles in pressed PBX 9501. The models of the dry blend have been labeled "DB" 
while those of pressed PBX 9501 have been labeled "PP". 

GMC calculations were performed on the PBX 9501 RVEs after discretizing each RVE into lOOx 100 subcells and 
assigning materials to subcells using the effective subcell approach. In order to validate the GMC predictions, FEM 
calculations were also performed on the RVEs after discretizing each RVE into 350 x 350 four-noded square elements. 
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The binary subcell approach was appHed to assign materials to elements for the FEM calculations. 

The particles in each RVE were assigned properties of HMX from Table |2] However, since particles occupy 
92% of the total volume in actual PBX 9501 while the sample microstructures could be filled only up to ~'86%, an 
intermediate homogenization step was required to determine the properties of the binder. To produce the desired 
92% volume of particles, a fine-particle filled binder containing 36% particles by volume, or "dirty" binder, was 
assumed. The effective elastic properties of the dirty binder were calculated using the differential effective medium 
approximation (Markov 2000). 

Table |4] shows the effective stiffness from FEM and GMC calculations for the models of PBX 9501 shown in 
Figure |5] For all microstructures, the values of Cf( and C|f predicted by GMC are less than 5% of the FEM values 
and less than 10% of the experimental values for PBX 9501 (shown in Table|2|. The FEM estimates increase with RVE 
size (varying from 150% to 450% of the experimental values), reflecting the dependence of predicted stiffnesses on 
microstructure, discretization and particle size distribution. The RVEs contain numerous particle to particle contacts, 
the number of which increases with increase in RVE size. These contacts lead to significant stress-bridging and hence 
relatively high values of stiffness as is reflected in the FEM predictions. However, stress-bridging is not incorporated 
accurately in the GMC approach leading to considerably lower values of effective stiffness. The issue of stress-bridging 
is further explored in the following section. 

The values of C|| predicted by GMC are around 0.5% of the FEM predictions. The reason for this large difference 
is that the effective shear stiffness predicted by GMC is simply the harmonic mean of the subcell shear moduli and 
only provides a lower bound on the shear stiffness. 

5 Stress-Bridging 

Comparisons of effective stiffness properties predicted by GMC with other numerical estimates have shown that GMC 
performs quite well for low modulus contrast materials with volume fractions below 60%. However, for high modulus 
contrast materials with high particle volume fractions, GMC usually predicts considerably lower effective stiffnesses 
than finite element calculations. In this section, GMC is applied to selected microstructures containing stress-bridging 
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and the predicted properties are compared to finite element estimates. The goal is to demonstrate that the effects of 
stress-bridging on effective properties are inaccurately described by GMC. 

5.1 Corner bridging : X-shaped microstructure 

In the RVE shown in Figure |6] the particles are square, arranged in the form of an 'X', and occupy a volume fraction 
of 25%. The particles transfer stress through corner contacts. The effective properties of the X-shaped microstructure 
shown in Figure|6]were calculated using the properties of HMX and five different binders with Young's moduli that 
range from 0.7 MPa to 7000 MPa, as shown in Table |5] 

Figure |2] shows the variation in the effective stiffness properties Cf[ and Cgf of the X-shaped microstructure with 
increasing Young's modulus contrast between the particles and the binder (Ep/Et,). These effective stiffness properties 
have been calculated using both finite elements (FEM) (256x256 elements) and GMC (64x64 subcells). The FEM 
and GMC estimates are in good agreement for Young's modulus contrasts of 200 or less. For higher Young's modulus 
contrasts, the effective stiffness properties predicted by GMC are much lower than those predicted by FEM. Note that 
the FEM estimates do not change significantly with increased discretization, implying that the solution has converged. 
The effect of corner singularities is also averaged out while calculating the effective properties using FEM. If it is 
assumed that the FEM estimates are close to the actual effective moduli of the RVE, the GMC estimates for high 
modulus contrasts are orders of magnitude lower than the actual effective moduli. Hence, GMC does not capture the 
stiffening effect of corner contacts accurately. 

Since the corner stress-bridging problem involves high stress concentrations that are not resolved well by finite 
elements, it is possible that the FEM calculations overestimate the effective properties of the X-shaped microstructure. 
Such corner singularities are minimized in the microstructures studied in the next section where the effect of stress- 
bridging along particles edges is studied. 

5.2 Edge bridging 

Figure[8]shows five RVEs (A through E) in which the degree of stress-bridging is increased progressively from corner 
bridging, to partial edge bridging, and finally to continuous stress-bridging across the RVE. In Figure H] the ' 1 ' 
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direction corresponds to the x-axis and the '2' direction corresponds to the y-axis. 

Model A contains a square particle that occupies 25% of the volume, is centered in the RVE, and does not have any 
stress-bridging. Model B contains three particles that contact along a diagonal of the RVE. In model C, particle contact 
is increased to produce a single line of stress-bridging in the a;-direction along the center of the RVE. Model D extends 
the line of contact in the x-direction to an area of contact in the a:-direction. In Model E, particle bridging across the 
RVE is extended to both directions. The material properties of the constituents of PBX 9501 at room temperature 
and low strain rate were used for the calculations (Table |2]l. GMC simulations of the RVEs were performed using 
lOOx 100 subcells while the validating finite element calculations were performed using approximately 10,000 eight- 
noded quadrilateral elements. Table|6]shows the effective stiffnesses of the five models obtained from GMC and finite 
element (FEM) calculations. 

As expected in model A, GMC and FEM predict nearly the same values of effective stiffness since there is no 
stress-bridging in the model (the effective stiffness is determined primarily by the volume fraction occupied by the 
square particle). However, FEM calculations for model B show that the diagonal stress-bridge in the model produces 
a higher stiffness than would occur if only the volume fraction occupied by the particles were considered in the 
calculation of effective stiffness. 

The GMC calculations for model B predict values of Cl^ and Cff that are lower than the FEM estimates by 
a factor of 18. This discrepancy implies that the diagonal stress-bridge in model B is not detected by the GMC 
calculations. The value of C|f from FEM is around 1,400 times that from GMC. This difference shows that, in the 
presence of stress-bridging, the shear stiffness can be considerably underestimated by GMC, even for low volume 
fraction composites. 

Model C has a continuous path through particles along the a;-axis (the ' 1 ' direction) and another continuous particle 
path along one diagonal. Intuitively, the stress-bridge path along the '1' direction is expected to primarily affect the 
normal components of stiffness (Cff, Cff, C|f) while particle contact along the diagonal is expected to affect the 
shear stiffness (C|f). These paths are shown by dashed lines (for normal stress-bridging) and by dotted lines (for 
shear stress-bridging) in Figure|9] Results for model C in Table|6]show that FEM predicts a considerable stiffening in 
the '1' direction while GMC does not appear to account for these stress-bridges. Since the shear stiffness from GMC 
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is simply a harmonic means of the subcell stiffnesses, C|f is not affected at all by geometry and only increases in 
proportion with the volume fraction of particles in the RVE. 

The estimates of C^f for model D and of Cff, Cff, C22 ^'^^ model E show that GMC can capture the effect 
of stress-bridging, provided there are continuous rows of particles with edge-to-edge contacts extending completely 
across the RVE. 

These studies of stress-bridging explain why GMC underestimates the effective modulus of the PBX 9501 models 
shown in Figure |5] In all these models, if 100x100 subcells are used to discretize the RVE, there are no rows or 
columns of subcells extending across the RVE that contain no binder Though corner contacts and other continuous 
stress paths exist in the PBX 9501 models, the effects of these stress-bridging paths are not incorporated into the GMC 
estimates of effective stiffness. The strain-compatible or shear-coupled method of cells (Williams & Aboudi 1999, Gan 
et al. 2000) approaches may be able to overcome some of these deficiencies of GMC. However, the computational 
efficiency of GMC is greatly reduced when these modifications are incorporated into GMC and hence the attractiveness 
of this micromechanics approach as an alternative to finite element analysis is also reduced. 

6 The recursive cell method 

The recursive cell method (RCM) (Banerjee 2002) is a real-space renormalization group (Wilson 1971, Wilson 1979) 
approach for calculating the effective elastic properties of composites that has been developed to address the short- 
comings of GMC while retaining high computational efficiency. A schematic of RCM is shown in Figure[TO] In RCM, 
as in GMC, the RVE is first discretized into a regular grid of subcells. For the first iteration of the recursive process, 
the subcells are assigned material properties based on the particle distribution in the RVE using the binary subcell ap- 
proach discussed earlier The subcells in the original grid are then grouped into blocks of nxn subcells. The effective 
elastic stiffness matrix of each of the blocks is calculated using a suitable homogenization approach such as GMC or 
FEM. Effective stiffnesses are assigned to each block, resulting in a new, coarser grid. This procedure is repeated until 
only one homogeneous block remains. The properties of this homogeneous block are the effective properties of the 
RVE. 
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Studies on the recursive cell method (Banerjee 2002) have shown that the method leads to an upper bound on the 
effective elastic properties if a FEM approach is used to homogenize blocks of subcells. As the number of elements 
used to discretized a block is increased, the value of the upper bound decreases and a more accurate estimate of the 
effective properties is obtained. GMC is an attractive alternative to the FEM approach for homogenization since less 
discretization is required to arrive at the same level of accuracy. 

In the previous section, GMC has been shown to not properly account for stress-bridging in the absence of con- 
tinuous stress-bridge paths across a RVE. However, error due to improper stress-bridging is reduced when GMC is 
used as the homogenizer in RCM because the probability of the existence of continuous stress-bridging paths across 
blocks of subcells is greater than that for the whole RVE. In addition, homogenization errors due to the overestimation 
or underestimation of stress-bridging in sections of the RVE are averaged out if the particle distribution is sufficiently 
random. 

A second source of error in GMC is the underestimation of the shear stiffness term Cgf . However, this error can be 
avoided while using RCM to determine the effective elastic properties of PBXs because, for macroscopically isotropic 
materials such as PBXs, relatively accurate estimates of the effective shear stiffness can be obtained from the effective 
normal stiffness terms (Banerjee 2002) and therefore direct estimates of Cgg are not required. On the other hand, if the 
composite is not macroscopically isotropic, a FEM homogenizer (Banerjee 2002) can be used to determine the value 
of C|f instead of GMC. 

The RCM technique has been applied to the four microstructures of the dry blend and pressed PBX 9501 shown in 
Figure lUa) and Elb), respectively. Each RVE was discretized into blocks of 256x256 square subcells of equal size. 
At each stage of recursion, blocks of 2 x2 subcells were homogenized using GMC. 

The values of Cjf for the four dry blend microstructures obtained from finite element (FEM) calculations, GMC 
calculations and RCM calculations are compared in Figure fTTT a). The RCM estimates of Cff for the four microstruc- 
tures vary from 90% to 150% of the FEM estimates. These RCM estimates are a considerable improvement over the 
GMC estimates shown as black bars in Figure fTTT a). Comparisons of Cjf for the four pressed PBX 9501 microstruc- 
tures (shown in Figure|5lb)) are shown in Figure fTTl b). For pressed PBX 9501, the RCM estimates vary between 84% 
and 180% of the finite element estimates. These RCM estimates are also a considerable improvement over the GMC 
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estimates of effective properties. RCM estimates of C22 and Cjf for the dry blend and pressed PBX 9501 have also 
been found to be in much better agreement with FEM results than the GMC estimates. 

As was expected, the estimated value of Cgg from RCM is quite low compared to both finite element estimates 
and experimental data. An improved estimate of Cgf can be obtained if the shear stiffness of each RCM block is 
calculated using finite elements (Banerjee 2002). The normal stiffnesses Cf(, Cjf, and C22 can be still be calculated 
using GMC, taking advantage of the absence of shear couphng. 

These results show that the RCM approach, in conjunction with a GMC homogenizer, can be used to arrive at 
reasonably accurate estimates of the effective properties of PBX materials. The RCM approach can therefore be 
used as an alternative to direct GMC calculations for high volume fraction, strong modulus contrast materials such as 
polymer bonded explosives. 

7 Summary and conclusions 

The generalized method of cells (GMC) has been found to accurately predict the effective elastic properties of com- 
posites containing square arrays of disks for volume fractions up to 0.60. However, for two-dimensional models of 
the polymer bonded explosive PBX 9501, estimates of effective elastic properties from GMC have been found to be 
considerably lower than both experimental values and estimates based on finite element (FEM) calculations. 

The lower values of normal stiffness predicted by GMC for PBX 9501 are due to inadequate incorporation of 
particle stress bridging into the approach. Model representative volume elements (RVEs) with comer and edge stress 
bridging show that corner bridging is ignored by GMC while edge stress bridging is incorporated only if continuous 
stress bridges exist along entire rows or colunms of subcells that traverse the length of the RVE. Low values of 
effective shear stiffness predicted by GMC can be attributed to the use of a harmonic mean of subcell shear stiffnesses 
to determine the effective shear stiffness of a RVE. The harmonic mean is a lower bound on the effective shear stiffness 
and is not apphcable for microstructures where there is significant interaction between particles. 

Improvements suggested to GMC that incorporate normal-shear coupling and strain compatibility across subcells 
have the potential to overcome some of these weaknesses of GMC. However, these improvements lead to much larger 
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systems of equations and a considerable increase in the computational cost of the method. The requirement of inverting 
a large matrix to obtain the effective properties makes the generalized method of cells very inefficient as the number 
of subcells increases. When materials such as PBX 9501 are modeled, the number of subcells needed to represent 
a random distribution of particles necessarily becomes large. In such situations, the generalized method of cells 
becomes inefficient and it may be preferable to perform finite element analyses to determine the effective properties. 
Thus, GMC does not appear to be an improvement over finite element analyses for high volume fraction, high modulus 
contrast particulate composites such as polymer bonded explosives. 

A computationally efficient alternative to both direct GMC and finite elements is the recursive cell method (RCM) 
with GMC being used to homogenize blocks of subcells. RCM estimates of normal stiffness terms for models of PBX 
9501 show considerable improvement compared to GMC estimates. The RCM estimates of shear stiffness can be 
improved if FEM is used, rather than GMC, to determine the effective shear stiffness of blocks of subcells. RCM, with 
a combination of GMC and FEM being used to homogenize blocks of subcells, has the potential of providing fast and 
accurate estimates of the effective properties of polymer bonded explosives. 
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Figure 1: RVE, subcells and notation used in GMC. 

Figure 2: Comparison of effective moduli of square arrays of disks from Greengard and Helsing (1998) (G&H) and 
GMC calculations. 

Figure 3: Manually generated microstructures containing ~ 90% circular particles by volume. 

Figure 4: Schematics of the application of the binary subceU approach and the effective subceU approach in GMC 
calculations. 

Figure 5: Microstructures containing circular particles based on the particle size distribution of the dry blend (DB) of 
PBX 9501 and of pressed (PP) PBX 9501. 

Figure 6: RVE used for comer stress-bridging model. 

Figure 7: Variation of effective stiffness with modulus contrast for 'X'-shaped microstructure. The Young's modulus 
contrast is the ratio of the Young's modulus of the particles to that of the binder. 

Figure 8: Progressive stress-bridging models A through E. 

Figure 9: Stress-bridging paths for Model C. 

Figure 10: Schematic of the recursive cell method. 

Figure 11: Comparisons of estimates of C^f for (a) models of the dry blend of PBX 9501 (b) models of pressed PBX 
9501. 
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Table 1: Component properties used by Greengard and Helsing (1998). 





Young's 


Poisson's 


Two-Dimensional 


Shear 




Modulus 


Ratio 


Bulk Modulus 


Modulus 




(MPa) 




(MPa) 


(MPa) 


Disks 


324 


0.20 


225 


135 


Binder 


2.7 


0.35 


3.3 


1 



Table 2: Experimentally determined elastic moduli and stiffness of PBX 9501 and its constituents (Wetzel 1999). Cij 
are components of the stiffness matrix. 



Material 


Young's 


Poisson's 


Cii = C22 


C12 


C'ee 




Modulus 


Ratio 










(MPa) 




(MPa) 


(MPa) 


(MPa) 


Particles 


15300 


0.32 


21894 


10303 


5795 


Binder 


0.7 


0.49 


11.97 


11.51 


0.235 


PBX 9501 


1013 


0.35 


1626 


875 


375 



Table 3: Effective stiffnesses of the six model microstructures from GMC and FEM calculations. 





Cfl (MPa) 


Cfi (MPa) 


Q« (MPa) 




FEM 


GMC 


FEM 


GMC 


FEM 


GMC 






Binary 


Effective 




Binary 


Effective 




Binary 


Effective 






Subcell 


Subcell 




Subcell 


Subcell 




Subcell 


Subcell 


Model 1 


177 


814 


479 


90 


119 


103 


11 


2.4 


2.3 


Model 2 


181 


807 


477 


86 


112 


103 


12 


2.3 


2.3 


Model 3 


186 


815 


193 


88 


108 


89 


15 


2.3 


2.2 


Model 4 


143 


116 


142 


114 


112 


124 


33 


2.4 


2.6 


Model 5 


237 


132 


323 


94 


100 


104 


38 


2.3 


2.5 


Model 6 


229 


132 


334 


76 


93 


100 


9 


2.1 


2.5 


Mean 


192 


471 


325 


91 


107 


104 


20 


2.3 


2.4 



Table 4: Effective stiffness of the model PBX 9501 microstructures from GMC and FEM calculations. 





Cfl (MPa) 


Cfi (MPa) 


Cfi (MPa) 


Q« (MPa) 




FEM 


GMC 


FEM 


GMC 


FEM 


GMC 


FEM 


GMC 


Model DBl 


2385 


152 


2094 


148 


633 


122 


792 


4.9 


Model DB2 


3618 


146 


1643 


144 


656 


122 


750 


4.9 


Model DB3 


3546 


149 


3385 


148 


1142 


125 


1317 


5.0 


Model DB4 


5274 


144 


5124 


146 


1712 


120 


1703 


4.8 


Model PPl 


3180 


180 


3570 


188 


989 


131 


1262 


4.8 


Model PP2 


3886 


170 


3683 


190 


1032 


132 


1278 


5.7 


Model PP3 


6302 


181 


6221 


181 


2043 


133 


2077 


6.6 


Model PP4 


7347 


182 


7587 


186 


2547 


129 


2542 


6.9 



20 



Table 5: The elastic properties of the components of the X shaped microstructure. Cij are components of the stiffness 
matrix. 





Young's 


Poisson's 


Cn 


Cl2 


Gee 




Modulus 


Ratio 










(MPa) 




(MPa) 


(MPa) 


(MPa) 


Particles 


15300 


0.32 


21894 


10303 


5795 


Binder a 


0.7 


0.49 


12 


11.5 


0.2 


Binder b 


7 


0.49 


120 


115 


2.4 


Binder c 


70 


0.49 


1198 


1151 


23.5 


Binder d 


700 


0.49 


11980 


11510 


235 


Binder e 


7000 


0.49 


119799 


115101 


2349 



Table 6: Effective properties of edge bridging models. 





Cfl (MPa) 


C^2 (MPa) 


Cfa' (MPa) 


Cl'i (MPa) 


FEM 


GMC 


FEM 


GMC 


FEM 


GMC 


FEM 


GMC 


Model A 


16 


16 


16 


16 


15 


15 


0.4 


0.3 


Model B 


336 


19 


343 


19 


337 


18 


537 


0.4 


Model C 


4095 


25 


889 


24 


1470 


23 


1093 


0.5 


Model D 


8992 


8540 


1361 


32 


523 


23 


1182 


0.6 


Model E 


10017 


9042 


10052 


9042 


2892 


2143 


1799 


0.9 
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Figure 1 
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Figure 2 
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Figure 3 
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Overlay Grid on RVE Assign Particle or Binder to Subcells 

(a) Binary subcell approach. 




Overlay Grid on RVE Assign Effective Properties to Subcells 

(b) Effective subcell approach. 

Figure 4 
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0.65x0.65mm2 0.94x0.94 mm^ l.Bxl.Bmm^ 1.33xl.33mm2 
DBl DB2 DBS DB4 

(a) Dry Blend of PBX 9501 




0.36x0.36nun2 0.42x0.42imn2 0.54x0.54 min^ 0.68x0.68 nun^ 
PPl PP2 PP3 PP4 

(b) Pressed PBX 9501 



Figure 5 
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Figure 7 
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Figure 10 
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